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•r 

Numerical Analysis - Approximating Functions 



Latest Modification: August 27, 2001 



Notation: Equation numbers are placed in parentheses, (Eq. 1), at the end of the sentence 
preceding the mathematical presentation of the equation. They should not be read as if they 
are part of the sentence. When reference is made to an equation, the word "equation" will be 
spelled out and not abbreviated. Mathematical symbols are not in all cases precise. Greek 
symbols will show as such on a Windows operating system, but not a Unix or Mac system. 



Series Approximation 

Given values of a function f(x) at a finite number of discrete points (n), 

x 0 ? x l 5 x 2' x 3' x n-l> x n 

f(x Q ), f( Xl ), f(x 2 ), f(x 3 ), fCx^), f(x n ) 

the goal is to construct an arbitrary, but simpler, function p(x) in the form of series that takes 
on the values of fj at the n points x i5 that is (Eq. 1), 



p(x) = f(x), for x Q < x < x n 
and possibly for x < Xq or x > x n 

that is to say (Eq. 2) 

POO = f(Xj) 5 for x 0 ,...,x n 

The function p(x) can take a variety of forms, such as 

• when p(x) is a polynomial, the process of representing f(x) by p(x) is called a 
polynomial approximation, 

• when p(x) is a finite trigonometric series, the process of representing f(x) by p(x) is 
called a trigonometric approximation, and 

• in like manner p(x) may be a series of: 

° exponential functions 

o Legendre polynomials 

o Bessel functions 

o etc. 



Approximating functions are the basis for interpolation schemes, and interpolation schemes 
are the basis for numerical quadrature relations. Thus approximating functions represent an 
important element in numerical methods. 
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Difference Tables 

Before we consider polynomial interpolation schemes, we should develop the concept of 
difference tables in the following document called Difference Tables . 



Interpolation 

Now that we have the concept of differences and their notation, we can develop the most 
elementary interpolation schemes known as polynomial interpolation found in the following 
document Polynomial Interpolation . 

A frequently encounter interpolation problem is one in which the data is not equally spaced 
in the independent variable. This leads us to consider the material in the following 
document on Non-Equidistant Interpolation . 



Least Square Data Analysis 

Finally, we should explore a different approach to approximating functions known as least 
squares as discussed in the following document Least Squares . 
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Difference Tables 



Latest Modification: September 28, 1999 



Notation: Equation numbers are placed in parentheses, (Eq. 1), at the end of the sentence 
preceding the mathematical presentation of the equation. They should not be read as if they 
are part of the sentence. When reference is made to an equation, the word "equation" will be 
spelled out and not abbreviated. Mathematical symbols are not in all cases precise. Greek 
symbols will show as such on a Windows operating system, but not a Unix or Mac system. 



Horizontal Difference Table 

Difference tables lend themselves to the use of spreadsheets as a means of solving many 
problems rather than resorting to writing computer code in either Fortran or C/C++. One 
type of difference table is known as a horizontal difference table in which differences are 
maintained on horizontal rows with the function as shown below for which successive 
differences through the fourth are denoted by Ajf^x), A 2 f i (x) 5 A 3 fj(x), and A 4 fj(x). 
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Error Propagation in Difference Tables 

Errors propagate in a horizontal difference table as shown below, where s is a small 
increment in the function f. 
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There are several points one should note about the propagation of the error. 

1. The error moves to successive rows with successive differences. 

2. The coefficients of s's are binomial coefficients with alternating signs. 

3. The algebraic sum of the e's in any difference column is zero. 

4. The first erroneous difference of any order is on the same horizontal line as the 
erroneous functional value. 
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Polynomial Interpolation 



Latest Modification: November 1, 2001 



Notation: Equation numbers are placed in parentheses, (Eq. 1), at the end of the sentence 
preceding the mathematical presentation of the equation. They should not be read as if they 
are part of the sentence. When reference is made to an equation, the word "equation" will be 
spelled out and not abbreviated. Mathematical symbols are not in all cases precise. Greek 
symbols will show as such on a Windows operating system, but not a Unix or Mac system. 



Polynomial Interpolation 

The justification for approximating with polynomials rest on a theorem by Weierstrass in 
1885, to which 

"...every continuous function in an interval (a 5 b) can be represented in that 
interval to any desired accuracy by a polynomial." 

One can also show that the following propositions are true for a polynomial 

1. nth differences of an nth-degree polynomial are constant for an arithmetic progression 
(equal interval) of the independent variable, 

2. if nth differences of a tabulated function for equal interval data are constant, the 
unknown function is a polynomial of degree n. 



Example 

Let f(x) = x 2 + 3x + 2. Forming differences in a horizontal difference table, we have the 
table shown below. 
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Since the function is indeed a 2nd degree polynomial the 2nd differences are constant. 



Polynomial Interpolation Retaining Second Differences 

If we restrict our consideration to the retention of second differences and ignore all 
higher-order differences, there are four useful polynomial approximations that can be 
derived for equal-interval interpolation. 

Newton's Forward Interpolation: Given values of a function f(x) at a finite number of 
discrete points (n), 



x 0> x l> x 2> x 3 5 — > x n-l> 



x 



n 



f ( x o) = f 0> f ( x l) = f l> f ( X 2> = f 2> f ( X 3) = f 3> -> f ( X n-l) = f n-l' f ( x n) = f n 
and letting u = [(x - Xj) / (x i+J - x^], then Newton's forward interpolation formula can be 
expressed using difference notation as 

p(x) = p(x 0 + uAx) 

or substituting the difference values of the function, we have 

p(x) = f Q + uAjfj + [u(u-l)/2!]A 2 f 2 + [u(u-l)(u.2)/3!]A 3 f 3 +...+ 
[u(u-l)(u-2)...(u-n+l)/n!]A n f n -H... 

which can be simplified by retaining only second order differences and substituting the 
functional values for the differences (Eq. 1) 

p(x) = f 0 + u(f 1 -f 0 ) + [u(u-l)/2 ](f 2 -2f!+f 0 ) 
Since Newton's forward interpolation formula uses values of the function at points ahead in 
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. the table, it can be used at the beginning of the tabulated values, but not at the end. 

Newton's Backward Interpolation: Given values of a function f(x) at a finite number of 
discrete points (n), and letting u = [(x - Xj) I (xj - x^)], then Newton's backward 

interpolation formula can be expressed using difference notation as 

p(x) = p(x Q + uAx) 

or substituting the difference values of the function, we have 

p(x) = f n + uA t f n + [u(u+l)/2!]A 2 f n + [u(u+l)(u+2)/3!]A 3 f n +...+ 
[u(u+l)(u+2)...(u+n-l)/n!]A n f n +... 

which can be simplified by retaining only second order differences and substituting the 
functional values for the differences (Eq. 2) 

p(x) = f 0 + u(f 0 -f 1 ) + [u(u+l)/2](f 0 -2f 1 +f 2 ) 

which is useful for interpolating near the end of the tabulated functions. 

Stirling's Central-Difference Interpolation: This interpolation formula is based on a 
diagonal difference table rather than a horizontal difference table. Without going in to the 
distinctions between the two types of difference tables, we can write down an interpolation 
formula including only second differences that is similar to Newton's interpolation formulas. 
Thus given values of a function f(x) at a finite number of discrete points (n), and letting u = 
[(x - Xj) / (x i+l - Xj)], Stirling's central-difference interpolation formula can be expressed 

using functional values as (Eq. 3) 

p(x) = f 0 + (u/2)(f r f 1 ) + [u 2 /2](f 1 -2f 0 + f 1 ) 

which also uses three successive points to approximate the function as does Newton. 

Bessel's Central-Difference Interpolation: This interpolation formula is also based on a 
diagonal difference table rather than a horizontal difference table. We shall also write down 
an interpolation formula including only second differences that is similar to Newton's and 
Stirling's interpolation formulas. Thus given values of a function f(x) at a finite number of 
discrete points (n), and letting u = [(x - Xj) / (x j+1 - Xj)], Bessel's central-difference 

interpolation formula can be expressed using functional values as (Eq. 4) 
p(x) = f 0 + u(f 1 -f 0 ) + [u(u-l)/4](f 2 -f 1 -f 0 + f 1 ) 

which uses four successive points rather than three to approximate the function. Note that in 
equations 1, 2, 3, and 4, there is a zero-order term, a first-order term, and a second-order 
term. The zero- and first-order terms constitute essential a linear interpolation, while the 
second-order term is a correction term to that linear interpolation. 



Example 
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4 Find f(x) = exp(-x), at the point x = 1 .7565, given the following values of the function. 
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Forming the value of u for the four equidistance interpolation formulas, we have for 
Newton's forward, Stirling's, and Bessel's interpolation formulas 

u = ( x - Xq ) / ( x x - x 0 ) = ( 1.7565 - 1.75 ) / ( 1.76 - 1.75 ) - 0.0065 / 0.01 = + 0.65 , 

Newton's Forward Interpolation: 

u ( u - 1 ) / 2 = 0.65 ( 0.65 - 1 ) / 2 = - 0.1 1375 , 
Stirling's Central-Difference Interpolation: 

u/2 = 0.65/2 = 0.325, 

u 2 /2 = (0.65 ) 2 /2 = 0.21125, 
Bessel's Central-Difference Interpolation: 

u ( u - 1 ) / 4 - 0.65 ( 0.65 - 1 ) / 4 = - 0.056875 . 
For Newton's backward formula, the value of u is given by 

u = ( x - Xj ) / ( x 1 - Xq ) = ( 1.7565 - 1.76 )/( 1.76 - 1.75 ) = - 0.0035 / 0.01 = - 0.35 , 

Newton's Backward Interpolation: 

u ( u + 1 ) / 2 = - 0.35 ( - 0.35 + 1 ) / 2 = - 0.1 1375 . 

The following table evaluates the three terms that appear in each of the four interpolation 
formula. Note that in order to apply Newton's backward formula, we must re-label the 
points so that i = -1, 0, 1 become i = -2, -1, 0. 
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It is evident that Bessel's interpolation formula has given us the most accurate result. This is 
inpart true because Bessel's formula is known to be more accurate when u lies between 0.25 
and 0.75. 



Non-Equidistant Interpolation 

A frequently encounter interpolation problem is one in which the data is not equally spaced 
in the independent variable. This leads us to consider Lagrangian interpolation in the 
following document on Non-Equidistant Interpolation . 
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Non-Equidistant Interpolation 

Lagrange's Formula 



Latest Modification: November 1, 2001 



Notation: Equation numbers are placed in parentheses, (Eq. 1), at the end of the sentence 
preceding the mathematical presentation of the equation. They should not be read as if they 
are part of the sentence. When reference is made to an equation, the word "equation" will be 
spelled out and not abbreviated. Mathematical symbols are not in all cases precise. Greek 
symbols will show as such on a Windows operating system, but not a Unix or Mac system. 



Non-Equidistant Data Interpolation 

The previous discussion of interpolation is predicated on equidistant intervals in the 
independent variable x. Sometimes it is inconvenient or impossible to obtain function 
values on an equidistant interval. In such cases, Lagrange's formula can be fashioned to only 
utilize such data as may be available. 

The general form of Lagrange's interpolation formula is given by (Eq. 1) 

P(x) = j=0 2 n { [ i=o nn ( x - xj ) ] / [ ( x - Xj ) i=0 • mt .Tl" ( Xj - Xj ) ] } 

where 

i=Q n n ( X - Xj ) = ( X - Xj ) ( X - Xj )...( x - x n ) 
i=0,i not j nn ( Xj - xj ) = ( Xj - x 0 )...( Xj - x M ) ( Xj - x i+1 )...( Xj - x n ) 
We should note that 

1 . lagrange's formula does not involve successive differences as does Newton's, 
Stirling's and Bessel's interpolation formulas. 

2. Lagrange's formula does involve n+1 successive pairs of variables (x n ,f n ). 

3. Although Lagrange's formula works for non-equidistant data, it is not restricted to 
such a limitation, i.e., it can be applied to equidistant data. 

4. If we know rates of change (differences) or rates of rates of change of the function, 
their neglect is to ignore information. 



Cubic Lagrangian Interpolation 

Suppose one wishes to construct a cubic polynomial through four successive points of a 
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. data set which we will label as (x^fj), (x 2 ,f 2 ), O^fjX ana (x 4 ,f 4 ). This can be accomplished 

from our general relation for n+1 points where we can construct the following relation (Eq. 
2) 



where 
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If one wishes, this expression can be written in the form (with a little algebra) as 

p(x) = ax 3 + bx 2 + cx + d 
However, from the computational point of view there is no real advantage in doing so. 



Example 

Find the cubic polynomial whose graph contains the four successive points (0,1), (1,2), 
(2,0), and (3,-2). Setting Xj = 0, fj = 1, x 2 = 1, f 2 = 2, x 3 = 2, f 3 = 0, and x 4 = 3, f 4 = -2, we 

can form the values of the Is 

1 14 = [ ( X - X 2 ) ( X - X 3 ) ( X - X 4 ) ] / [ ( Xj - X 2 ) ( Xj - x 3 ) ( Xj - x 4 ) ] 

=[(x-l)(x-2)(x-3)]/[(0-l)(0-2)(0-3)] 
=[(x-l)(x-2)(x-3)]/-6 

1 2 ,4 = [ ( x - Xj ) ( x - x 3 ) ( x - x 4 ) ] / [ ( x 2 - x 1 ) ( x 2 - x 3 ) ( x 2 - x 4 ) ] 
=[(x-0)(x-2)(x-3)]/[(l-0)(l-2)(l-3)] 
=[x(x-2)(x-3)]/2 

1 3 4 = [ ( X - Xj ) ( x - x 2 ) ( x - x 4 ) ] / [ ( x 3 - x, ) ( x 3 - x 2 ) ( x 3 - x 4 ) ] 

=[(x-0)(x-l)(x-3)]/[(2-0)(2-l)(2-3)] 
=[x(x-l)(x-3)]/-2 

1 4 4 = [ ( X - Xj ) ( x - x 2 ) ( x - x 3 ) ] / [ ( x 4 - x, ) ( x 4 - x 2 ) ( x 4 - x 3 ) ] 

=[(x-0)(x-l)(x-2)]/[(3-0)(3-l)(3-2)] 
=[x(x-l)(x-2)]/6. 

Inserting the values for the Is into Equation 2, we have 

f ( x ) = 1 l,4 f l +1 2,4 f 2 + 1 3,4 f 3 +1 4,4 f 4 
={[(x-l)(x-2)(x-3)]/-6}l+{[x(x-2)(x-3)]/2}2+{[x(x-l)(x- 
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. . 3)]/-2}0+{[x(x-l)(x-2)]/6}-2 

= (- 1 / 6 )[(x-l)(x-2)(x-3)] + [x(x-2)(x-3)] + (-V 3 )[x(x-l)(x-2)]. 

Multiplying out the terms and collecting yields the desired polynomial 

f(x) = V 2 x 3 - 3x 2 + 7 / 2 x + 1 . 
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Least Squares 



Latest Modification: November 14, 2001 



Notation: Equation numbers are placed in parentheses, (Eq. 1), at the end of the sentence 
preceding the mathematical presentation of the equation. They should not be read as if they 
are part of the sentence. When reference is made to an equation, the word "equation" will be 
spelled out and not abbreviated. Mathematical symbols are not in all cases precise. Greek 
symbols will show as such on a Windows operating system, but not a Unix or Mac system. 



Defining a "Close Fitting 1 ' Approximating Function 

Earlier we had discussed approximating functions ( Approximating Functions ). In that 
discussion the approximating function contained all of the points of a given discrete 
function, i.e., the function was constructed such that it agreed precisely with the discrete 
points used in its construction. Suppose now we wish to construct approximating functions 
whose graphs need not necessarily contain the points of given discrete function, but which 
"fits the function closely." Such is the method known as least squares. 

In order to understand this idea of "fits the function closely," let us consider a simple 
example. Given the points ( 0, -1 ), ( 1, 1 ), ( 2, 3 ), and ( 3, 4 ), let us find the straight line 
f(x) = aj + a2 x 5 which in some sense "fits the data closely." 



Figure 1 shows a graph of 
the four points that we wish 
to fit. One means of 
determining what is a close 
fit, is to minimize the 
vertical distances between 
the given points and the 
constructed line. Let is 
calculate the "errors" or 
"how far the line misses the 
given points" as follows 
(Eq. 1) 

8 1 = a ^ + &2 x i " fj > 
j>2 a | a^ x^ ~ , 

s 3 = a l + ^2 x 3 ' *3 ' 
e 4 = a l + *2 X 4 " f 4 * 

where our object is to minimize the total error. We could do that by considering the total 
error E. Let us first consider just the sum of the individual errors as given by (Eq. 2) 
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However, a simple sum will not suffice as the following example illustrates. Suppose the 
values of the errors are as shown 

E = ( 0.1 ) + ( -0.1 ) + ( 0.2 ) + ( -0.2 ) = 0 . 

Although the sum gives us a minimum, it is meaningless since clearly the total error is not 
zero. Another possible means of minimizing is to take the absolute values of the errors, but 
this means that the function E is not always differentiable. To preserve the differentiability 
of E and yet avoid the sign problem, we can minimize the squares of the errors as in (Eq. 3) 

E = e 1 2 + e 2 2 + e 3 2 + e 4 2 . 

Equation 3 is the desirable definition of E since (i) it is always differentiable, (ii) individual 
terms are non-negative, and (iii) because of the parabolic structure E always has a unique 
minimum. 



Minimizing the Total Error 

Let us construct our total error function for the given example (Eq. 4) 

E = [ ai -f l] 2 + [(a 1 + a 2 )-l] 2 + [(a 1 -f2a 2 )-3] 2 -f [(a 1 +3a 2 )-4] 2 . 

In order to minimize, we want to set the partial derivative of E with respect to each of the 
parameters in the linear fit, aj and a2> as shown (Eq. 5) 

5E / 5aj = 0 , 8E / bs^ = 0 . 

which yields the two linear algebraic equations 

4 aj + 6 a2 = 7 , 
6aj + 14^ = 19, 

whose solution is 

aj =-0.8 , 
82=1.7. 
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,Jm\s the approximating 
v function that "fits the data 
closely" is (Eq. 6) 

f(x) = a t + a 2 x = -0.8 + 1.7 



and the sense in which it fits 
closely is that it minmizes 
the squares of the vertical 
distances between the 
approximating function and 
the discrete points. 




imalipn/leastSquares.htm 1 



Figure 2 shows a graph of 
the four points that we wish 

to fit along with approximating function that minimizes the squares of the vertical distances 
between the given points and the constructed line. 



Are there other ways to define "fits closely?" Yes, for example, minimize the squares of the 
perpendicular distances rather than the vertical distances. However, this is a considerably 
harder problem to solve. Therefore, we will restrict ourselves to the standard definition of 
the least squares technique and minimize the vertical distances. 



Periodic or Quasi-Periodic Approximating Functions in Least Squares. 

Let us consider a second example, which reveals more of the subtle aspects of the method 
of least squares. 

the general discription of the method of least squares is as follows. Given the set of discrete 
points ( Xj, fj ), ( x 2 , f 2 ), ( x 3 , f 3 ), ... ( x n , f n ), let g(x) be a continuous function chosen 

such as to minimize the squares of the vertical distances, that is, minimize the total error E 
(Eq.7) 

E = [g(x 1 )-f 1 ] 2 + [g(x 2 )-f 2 )] 2 + [g(x3-f 3 ] 2 + ... + [g(x n -f n ] 2 . 

assuming that g(x) contains m independent parameters aj, e^, a m . In order to minimize E 

with respect to each of the m independent parameters, we solve the system of algebraic 
equations given by (Eq. 8) 

5E / 5aj = 0 , 8E / hd^ = 0 , 5E / 5a 3 - 0 , 5E / 5a m = 0. 
to yield the desired values of the m independent parameters a l5 e^, a m . 
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> >4ow let us consider the 
' * oscilloscope trace shown in 
* Figure 3. For our 

oscilloscope trace, what is 
the functional form for g(x) 
that "fits closely" the data? 
Some possibilities are (Eq. 
9) 

linear: g(x) = aj + x, 

quadratic: g(x) = aj + x + 
2 

a 3 X > 

cubic: g(x) = aj + ?^ x + a 3 
x^ + a 4 x . 
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However, the choice of the approximating function g(x) to be used must depend on some 
knowledge of the behavior of the phenomenon under observation. For the oscilloscope trace 
in Figure 3, the behavior is clearly periodic and a function like (Eq. 10) 

g(x) = a t sin(x) + cos(x) . 

would be a reasonable place to begin. Let us assume that the five points ( 0, -0.9 ) 5 ( 7i/4, 1 .5 
), {nil, 3.1 ), 37r/4, 3.0 ), and (rc, 1.1 ) are on the graph in Figure 3, so that the total error is 

E = [ aj sin(0) + ^ cos(0) + 0.9 ] 2 + [ aj sin(7t/4) + ^ cos(tt/4) - 1 .5 ] 2 + [ ^ sin(7c/2) + ^ 
cos(tt/2) - 3.1 ] 2 + [ aj sin(37r/4) + ^ cos(37t/4) - 3.0 ] 2 + [ a 2 sin(7c) + cos(tt) - 1.1 ] 2 . 



Next minimizing with 
respect to the two 
independent parameters aj 

and &2 yields two algebraic 

expressions that can 
immediately be solved for aj 

and ^ or 
aj = 3.14 and a2 = -1.02, 

so that (Eq. 11) 

g(x) = 3.14 sin(x)- 1.02 
cos(x) , 




which is shown in Figure 4 as the magenta curve. 
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